/*
 * =====================================================================================
 *
 *       Filename:  thin.cpp
 *
 *    Description:  Defines function that returns thinned version of binary image
 *                  Thinning algorithm by Guo & Hall
 *
 *        Version:  1.0
 *        Created:  19/07/12 11:00:12
 *       Revision:  none
 *       Compiler:  gcc
 *
 *         Author:  Moos Hueting (mh), mooshueting@gmail.com
 *   Organization:  
 *
 * =====================================================================================
 */
#include <iostream>
#include <cv.h>
#include <highgui.h>
#include "thin.h"

using namespace std;
using namespace cv;

// thinning algorithm by Guo and Hall
Mat thin(Mat src) {

    // keep thinning src until there is no change
    while (thinOneIt(&src))
        ;

    return src;

}

int thinOneIt(Mat* src) {
    Size s = src->size();

    Mat mask(s.height, s.width, CV_8UC1, Scalar(1));
    int unchanged = 0; // number of unchanged pixels

    // set multithreading using OpenMP
//#pragma omp parallel for shared(mask) schedule(dynamic, 1)
    // first subiteration
    for (int i=1; i < s.height - 1; i++) {
        for (int j=1; j < s.width - 1; j++) {
            // get field index
            int c = s.width*i + j;
            
            // if pixel is set, check conditions
            // if all met, delete from matrix
            if (src->at<unsigned char>(c) != 0) {
                // condition G1
                int xh = 0;
                if (src->at<unsigned char>(c + 1) == 0) {
                    if (src->at<unsigned char>((c - s.width) + 1) != 0 || src->at<unsigned char>(c - s.width) != 0)
                        ++xh;
                }
                if (src->at<unsigned char>(c - s.width) == 0) {
                    if (src->at<unsigned char>(c - (s.width + 1)) != 0 || src->at<unsigned char>(c - 1) != 0)
                        ++xh;
                }
                if (src->at<unsigned char>(c - 1) == 0) {
                    if (src->at<unsigned char>(c + (s.width - 1)) != 0 || src->at<unsigned char>(c + s.width) != 0)
                        ++xh;
                }
                if (src->at<unsigned char>(c + s.width) == 0) {
                    if (src->at<unsigned char>(c + s.width + 1) != 0 || src->at<unsigned char>(c + 1) != 0)
                        ++xh;
                }
                
                if (xh != 1)
                    continue;

                // condition G2
                int n1 = 0;
                int n2 = 0;

                if (src->at<unsigned char>(c + 1) != 0 || src->at<unsigned char>((c - s.width) + 1) != 0)
                    ++n1;
                if (src->at<unsigned char>(c - s.width) != 0 || src->at<unsigned char>(c - (s.width + 1)) != 0)
                    ++n1;
                if (src->at<unsigned char>(c - 1) != 0 || src->at<unsigned char>(c + (s.width - 1)) != 0)
                    ++n1;
                if (src->at<unsigned char>(c + s.width) != 0 || src->at<unsigned char>(c + s.width + 1) != 0)
                    ++n1;

                if (src->at<unsigned char>((c - s.width) + 1) != 0 || src->at<unsigned char>(c - s.width) != 0)
                    ++n2;
                if (src->at<unsigned char>(c - (s.width + 1)) != 0 || src->at<unsigned char>(c - 1) != 0)
                    ++n2;
                if (src->at<unsigned char>(c + (s.width - 1)) != 0 || src->at<unsigned char>(c + (s.width - 1)) != 0)
                    ++n2;
                if (src->at<unsigned char>(c + s.width + 1) != 0 || src->at<unsigned char>(c + s.width) != 0)
                    ++n2;

                if (n1 < n2) {
                    if (!(2 <= n1 && n1 <= 3))
                        continue;
                } else {
                    if (!(2 <= n2 && n2 <= 3))
                        continue;
                }

                // condition G3
                if ((src->at<unsigned char>((c - s.width) + 1) || src->at<unsigned char>(c - s.width) || !(src->at<unsigned char>(c + s.width + 1))) && src->at<unsigned char>(c + 1))
                    continue;

                // still here? then all conditions are met and we can delete the pixel
                mask.at<unsigned char>(c) = 0;
            }
        }
    }

    // apply mask
    multiply(*src, mask, *src);
    unchanged += countNonZero(mask);
    mask = Mat(s.height, s.width, CV_8UC1, Scalar(1));

    // second subiteration, only odd fields
//#pragma omp parallel for schedule(dynamic, 1)
    for (int i=1; i < s.height - 1; i++) {
        for (int j=1; j < s.width - 1; j++) {
            // get field index
            int c = s.width*i + j;

            // if pixel is set, check conditions
            // if all met, delete from matrix
            if (src->at<unsigned char>(c) != 0) {
                // condition G1
                int xh = 0;
                if (src->at<unsigned char>(c + 1) == 0) {
                    if (src->at<unsigned char>((c - s.width) + 1) != 0 || src->at<unsigned char>(c - s.width) != 0)
                        ++xh;
                }
                if (src->at<unsigned char>(c - s.width) == 0) {
                    if (src->at<unsigned char>(c - (s.width + 1)) != 0 || src->at<unsigned char>(c - 1) != 0)
                        ++xh;
                }
                if (src->at<unsigned char>(c - 1) == 0) {
                    if (src->at<unsigned char>(c + (s.width - 1)) != 0 || src->at<unsigned char>(c + s.width) != 0)
                        ++xh;
                }
                if (src->at<unsigned char>(c + s.width) == 0) {
                    if (src->at<unsigned char>(c + s.width + 1) != 0 || src->at<unsigned char>(c + 1) != 0)
                        ++xh;
                }
                
                if (xh != 1)
                    continue;

                // condition G2
                int n1 = 0;
                int n2 = 0;

                if (src->at<unsigned char>(c + 1) != 0 || src->at<unsigned char>((c - s.width) + 1) != 0)
                    ++n1;
                if (src->at<unsigned char>(c - s.width) != 0 || src->at<unsigned char>(c - (s.width + 1)) != 0)
                    ++n1;
                if (src->at<unsigned char>(c - 1) != 0 || src->at<unsigned char>(c + (s.width - 1)) != 0)
                    ++n1;
                if (src->at<unsigned char>(c + s.width) != 0 || src->at<unsigned char>(c + s.width + 1) != 0)
                    ++n1;

                if (src->at<unsigned char>((c - s.width) + 1) != 0 || src->at<unsigned char>(c - s.width) != 0)
                    ++n2;
                if (src->at<unsigned char>(c - (s.width + 1)) != 0 || src->at<unsigned char>(c - 1) != 0)
                    ++n2;
                if (src->at<unsigned char>(c + (s.width - 1)) != 0 || src->at<unsigned char>(c + (s.width - 1)) != 0)
                    ++n2;
                if (src->at<unsigned char>(c + s.width + 1) != 0 || src->at<unsigned char>(c + s.width) != 0)
                    ++n2;

                if (n1 < n2) {
                    if (!(2 <= n1 && n1 <= 3))
                        continue;
                } else {
                    if (!(2 <= n2 && n2 <= 3))
                        continue;
                }

                // condition G3'
                if ((src->at<unsigned char>(c + (s.width - 1)) || src->at<unsigned char>(c + s.width) || !(src->at<unsigned char>(c - (s.width + 1)))) && src->at<unsigned char>(c - 1))
                    continue;

                // still here? then all conditions are met and we can delete the pixel
                mask.at<unsigned char>(c) = 0;
            }
        }
    }

    unchanged += countNonZero(mask);
    multiply(*src, mask, *src);
    return 2*s.width*s.height - unchanged;

}
